I want to evaluate the function \[
f(x, y) = \frac{1}{|0.5 - x^4 - y^4| + 0.1}
\] which is not particularly well-behaved in that it has a discontinuity due to the absolute value, and it does not disappear at the boundaries. Evaluating the function through classical sparse grids is insufficient because ot requires too many points to capture the ridge. Adaptive sparse grids are better suited for this problem. This write-up demonstrates the potential of adaptive sparse grids.
Brute Force
First, I directly evaluate the function over the unit square, \([0, 1]^2\) using \(100^2\) points.
library('plotly')
data.classic <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/CSG.csv', header = FALSE)
fig1 <- plot_ly(data.classic, x = ~V1, y = ~V2, z = ~V3, marker = list(color = ~V3, colorscale = 'Portland', showscale = TRUE, size = 3))
fig1 <- fig1 %>% layout(title = 'Direct Evaluation')
fig1 <- fig1 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'f(x, y)')))
fig1
Classical Sparse Grids
Now I evaluate the function on sparse grids of level 14. In two dimensions, that means almost 280,000 grid points.
fig2 <- plot_ly(data.classic, x = ~V1, y = ~V2, z = ~V4, marker = list(color = ~V4, colorscale = 'Portland', showscale = TRUE, size = 3))
fig2 <- fig2 %>% layout(title = 'Classical Sparse Grid')
fig2 <- fig2 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'f(x, y)')))
fig2
The error between brute force and sparse grids, calculated as the \(L_2\) norm is 0.0889.
\[
L_2 = \frac{1}{N} \bigg(\sum_{i = 1}^N |f(x_i) - u(x_i)|^2\bigg)^\frac{1}{2}
\]
Adaptive Sparse Grids
Finally, I evaluate the function with adaptive sparse grids with level 5. In two dimensions, the algo starts with 257 bounds (boundary included) and iterates until the error is less than 0.1 which happens with 6571 grid points and the error is 0.0897.
data.adaptive <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/ASG.csv', header = FALSE)
fig3 <- plot_ly(data.adaptive, x = ~V1, y = ~V2, z = ~V4, marker = list(color = ~V4, colorscale = 'Portland', showscale = TRUE, size = 3))
fig3 <- fig3 %>% layout(title = 'Adaptive Sparse Grids')
fig3 <- fig3 %>% layout(scene = list(xaxis = list(title = 'x'),
yaxis = list(title = 'y'),
zaxis = list(title = 'f(x, y)')))
fig3
Life Cycle Problem
Basic
data.life <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/Life.csv', header = FALSE)
fig4 <- plot_ly(data.life, x = ~V1, y = ~V2, z = ~V3, marker = list(color = ~V3, colorscale = 'Portland', showscale = TRUE, size = 3))
fig4 <- fig4 %>% layout(title = 'Direct')
fig4 <- fig4 %>% layout(scene = list(xaxis = list(title = 'e'),
yaxis = list(title = 'x'),
zaxis = list(title = 'V(e, x)')))
fig4
fig5 <- plot_ly(data.life, x = ~V1, y = ~V2, z = ~V4, marker = list(color = ~V4, colorscale = 'Portland', showscale = TRUE, size = 3))
fig5 <- fig5 %>% layout(title = 'Adaptive')
fig5 <- fig5 %>% layout(scene = list(xaxis = list(title = 'e'),
yaxis = list(title = 'x'),
zaxis = list(title = 'V(e, x)')))
fig5
With Nelder-Mead
data.lifeB <- read.csv('/home/ahyan/Dropbox/Housing Market and Leverage Cycle/Code/src/Sparse/LifeB.csv', header = FALSE)
fig6 <- plot_ly(data.lifeB, x = ~V1, y = ~V2, z = ~V3, marker = list(color = ~V3, colorscale = 'Portland', showscale = TRUE, size = 3))
fig6 <- fig6 %>% layout(title = 'Direct')
fig6 <- fig6 %>% layout(scene = list(xaxis = list(title = 'e'),
yaxis = list(title = 'x'),
zaxis = list(title = 'V(e, x)')))
fig6
LS0tCnRpdGxlOiAiU3BhcnNlIEdyaWRzIgphdXRob3I6IEFoeWFuIFBhbmp3YW5pCmRhdGU6IEF1Z3VzdCAxMXRoLCAyMDIwCm91dHB1dDogaHRtbF9ub3RlYm9vawotLS0KCkkgd2FudCB0byBldmFsdWF0ZSB0aGUgZnVuY3Rpb24KJCQKZih4LCB5KSA9IFxmcmFjezF9e3wwLjUgLSB4XjQgLSB5XjR8ICsgMC4xfQokJAp3aGljaCBpcyBub3QgcGFydGljdWxhcmx5IF93ZWxsLWJlaGF2ZWRfIGluIHRoYXQgaXQgaGFzIGEgZGlzY29udGludWl0eSBkdWUgdG8gdGhlIGFic29sdXRlIHZhbHVlLCBhbmQKaXQgZG9lcyBub3QgZGlzYXBwZWFyIGF0IHRoZSBib3VuZGFyaWVzLiBFdmFsdWF0aW5nIHRoZSBmdW5jdGlvbiB0aHJvdWdoIGNsYXNzaWNhbCBzcGFyc2UgZ3JpZHMgaXMgaW5zdWZmaWNpZW50CmJlY2F1c2Ugb3QgcmVxdWlyZXMgdG9vIG1hbnkgcG9pbnRzIHRvIGNhcHR1cmUgdGhlIHJpZGdlLiBBZGFwdGl2ZSBzcGFyc2UgZ3JpZHMgYXJlIGJldHRlciBzdWl0ZWQgZm9yIHRoaXMKcHJvYmxlbS4gVGhpcyB3cml0ZS11cCBkZW1vbnN0cmF0ZXMgdGhlIHBvdGVudGlhbCBvZiBhZGFwdGl2ZSBzcGFyc2UgZ3JpZHMuCgojIyBCcnV0ZSBGb3JjZQoKRmlyc3QsIEkgZGlyZWN0bHkgZXZhbHVhdGUgdGhlIGZ1bmN0aW9uIG92ZXIgdGhlIHVuaXQgc3F1YXJlLCAkWzAsIDFdXjIkIHVzaW5nICQxMDBeMiQgcG9pbnRzLgoKYGBge3Igd2FybmluZz1GQUxTRSwgbWVzc2FnZT1GQUxTRSwgZmlnLmFsaWduPSdjZW50ZXInfQoKbGlicmFyeSgncGxvdGx5JykKCmRhdGEuY2xhc3NpYyA8LSByZWFkLmNzdignL2hvbWUvYWh5YW4vRHJvcGJveC9Ib3VzaW5nIE1hcmtldCBhbmQgTGV2ZXJhZ2UgQ3ljbGUvQ29kZS9zcmMvU3BhcnNlL0NTRy5jc3YnLCBoZWFkZXIgPSBGQUxTRSkKCmZpZzEgPC0gcGxvdF9seShkYXRhLmNsYXNzaWMsIHggPSB+VjEsIHkgPSB+VjIsIHogPSB+VjMsIG1hcmtlciA9IGxpc3QoY29sb3IgPSB+VjMsIGNvbG9yc2NhbGUgPSAnUG9ydGxhbmQnLCBzaG93c2NhbGUgPSBUUlVFLCBzaXplID0gMykpCgpmaWcxIDwtIGZpZzEgJT4lIGxheW91dCh0aXRsZSA9ICdEaXJlY3QgRXZhbHVhdGlvbicpCgpmaWcxIDwtIGZpZzEgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ3gnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd5JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnZih4LCB5KScpKSkKCmZpZzEKCmBgYAoKIyMgQ2xhc3NpY2FsIFNwYXJzZSBHcmlkcwoKTm93IEkgZXZhbHVhdGUgdGhlIGZ1bmN0aW9uIG9uIHNwYXJzZSBncmlkcyBvZiBsZXZlbCAxNC4gSW4gdHdvIGRpbWVuc2lvbnMsIHRoYXQgbWVhbnMgYWxtb3N0IDI4MCwwMDAgZ3JpZCBwb2ludHMuCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CmZpZzIgPC0gcGxvdF9seShkYXRhLmNsYXNzaWMsIHggPSB+VjEsIHkgPSB+VjIsIHogPSB+VjQsIG1hcmtlciA9IGxpc3QoY29sb3IgPSB+VjQsIGNvbG9yc2NhbGUgPSAnUG9ydGxhbmQnLCBzaG93c2NhbGUgPSBUUlVFLCBzaXplID0gMykpCgpmaWcyIDwtIGZpZzIgJT4lIGxheW91dCh0aXRsZSA9ICdDbGFzc2ljYWwgU3BhcnNlIEdyaWQnKQoKZmlnMiA8LSBmaWcyICU+JSBsYXlvdXQoc2NlbmUgPSBsaXN0KHhheGlzID0gbGlzdCh0aXRsZSA9ICd4JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAneScpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgemF4aXMgPSBsaXN0KHRpdGxlID0gJ2YoeCwgeSknKSkpCgpmaWcyCmBgYAoKVGhlIGVycm9yIGJldHdlZW4gYnJ1dGUgZm9yY2UgYW5kIHNwYXJzZSBncmlkcywgY2FsY3VsYXRlZCBhcyB0aGUgJExfMiQgbm9ybSBpcyAwLjA4ODkuCgokJApMXzIgPSBcZnJhY3sxfXtOfSBcYmlnZyhcc3VtX3tpID0gMX1eTiB8Zih4X2kpIC0gdSh4X2kpfF4yXGJpZ2cpXlxmcmFjezF9ezJ9CiQkCgojIyBBZGFwdGl2ZSBTcGFyc2UgR3JpZHMKCkZpbmFsbHksIEkgZXZhbHVhdGUgdGhlIGZ1bmN0aW9uIHdpdGggYWRhcHRpdmUgc3BhcnNlIGdyaWRzIHdpdGggbGV2ZWwgNS4gSW4gdHdvIGRpbWVuc2lvbnMsIHRoZSBhbGdvIHN0YXJ0cwp3aXRoIDI1NyBib3VuZHMgKGJvdW5kYXJ5IGluY2x1ZGVkKSBhbmQgaXRlcmF0ZXMgdW50aWwgdGhlIGVycm9yIGlzIGxlc3MgdGhhbiAwLjEgd2hpY2ggaGFwcGVucyB3aXRoIDY1NzEKZ3JpZCBwb2ludHMgYW5kIHRoZSBlcnJvciBpcyAwLjA4OTcuCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CmRhdGEuYWRhcHRpdmUgPC0gcmVhZC5jc3YoJy9ob21lL2FoeWFuL0Ryb3Bib3gvSG91c2luZyBNYXJrZXQgYW5kIExldmVyYWdlIEN5Y2xlL0NvZGUvc3JjL1NwYXJzZS9BU0cuY3N2JywgaGVhZGVyID0gRkFMU0UpCgpmaWczIDwtIHBsb3RfbHkoZGF0YS5hZGFwdGl2ZSwgeCA9IH5WMSwgeSA9IH5WMiwgeiA9IH5WNCwgbWFya2VyID0gbGlzdChjb2xvciA9IH5WNCwgY29sb3JzY2FsZSA9ICdQb3J0bGFuZCcsIHNob3dzY2FsZSA9IFRSVUUsIHNpemUgPSAzKSkKCmZpZzMgPC0gZmlnMyAlPiUgbGF5b3V0KHRpdGxlID0gJ0FkYXB0aXZlIFNwYXJzZSBHcmlkcycpCgpmaWczIDwtIGZpZzMgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ3gnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd5JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnZih4LCB5KScpKSkKCmZpZzMKYGBgCgoKIyMgTGlmZSBDeWNsZSBQcm9ibGVtCgpCYXNpYwpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CgpkYXRhLmxpZmUgPC0gcmVhZC5jc3YoJy9ob21lL2FoeWFuL0Ryb3Bib3gvSG91c2luZyBNYXJrZXQgYW5kIExldmVyYWdlIEN5Y2xlL0NvZGUvc3JjL1NwYXJzZS9MaWZlLmNzdicsIGhlYWRlciA9IEZBTFNFKQoKZmlnNCA8LSBwbG90X2x5KGRhdGEubGlmZSwgeCA9IH5WMSwgeSA9IH5WMiwgeiA9IH5WMywgbWFya2VyID0gbGlzdChjb2xvciA9IH5WMywgY29sb3JzY2FsZSA9ICdQb3J0bGFuZCcsIHNob3dzY2FsZSA9IFRSVUUsIHNpemUgPSAzKSkKCmZpZzQgPC0gZmlnNCAlPiUgbGF5b3V0KHRpdGxlID0gJ0RpcmVjdCcpCgpmaWc0IDwtIGZpZzQgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ2UnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd4JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnVihlLCB4KScpKSkKCmZpZzQKCmBgYAoKCgpgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CgpmaWc1IDwtIHBsb3RfbHkoZGF0YS5saWZlLCB4ID0gflYxLCB5ID0gflYyLCB6ID0gflY0LCBtYXJrZXIgPSBsaXN0KGNvbG9yID0gflY0LCBjb2xvcnNjYWxlID0gJ1BvcnRsYW5kJywgc2hvd3NjYWxlID0gVFJVRSwgc2l6ZSA9IDMpKQoKZmlnNSA8LSBmaWc1ICU+JSBsYXlvdXQodGl0bGUgPSAnQWRhcHRpdmUnKQoKZmlnNSA8LSBmaWc1ICU+JSBsYXlvdXQoc2NlbmUgPSBsaXN0KHhheGlzID0gbGlzdCh0aXRsZSA9ICdlJyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB5YXhpcyA9IGxpc3QodGl0bGUgPSAneCcpLAogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgemF4aXMgPSBsaXN0KHRpdGxlID0gJ1YoZSwgeCknKSkpCgpmaWc1CmBgYAoKV2l0aCBOZWxkZXItTWVhZApgYGB7ciB3YXJuaW5nPUZBTFNFLCBtZXNzYWdlPUZBTFNFLCBmaWcuYWxpZ249J2NlbnRlcid9CgpkYXRhLmxpZmVCIDwtIHJlYWQuY3N2KCcvaG9tZS9haHlhbi9Ecm9wYm94L0hvdXNpbmcgTWFya2V0IGFuZCBMZXZlcmFnZSBDeWNsZS9Db2RlL3NyYy9TcGFyc2UvTGlmZUIuY3N2JywgaGVhZGVyID0gRkFMU0UpCgpmaWc2IDwtIHBsb3RfbHkoZGF0YS5saWZlQiwgeCA9IH5WMSwgeSA9IH5WMiwgeiA9IH5WMywgbWFya2VyID0gbGlzdChjb2xvciA9IH5WMywgY29sb3JzY2FsZSA9ICdQb3J0bGFuZCcsIHNob3dzY2FsZSA9IFRSVUUsIHNpemUgPSAzKSkKCmZpZzYgPC0gZmlnNiAlPiUgbGF5b3V0KHRpdGxlID0gJ0RpcmVjdCcpCgpmaWc2IDwtIGZpZzYgJT4lIGxheW91dChzY2VuZSA9IGxpc3QoeGF4aXMgPSBsaXN0KHRpdGxlID0gJ2UnKSwKICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgIHlheGlzID0gbGlzdCh0aXRsZSA9ICd4JyksCiAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICB6YXhpcyA9IGxpc3QodGl0bGUgPSAnVihlLCB4KScpKSkKCmZpZzYKCmBgYAo=